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ABSTRACT 

A new computational scheme is developed to study gas accretion from a 
circumbinary disk. The scheme decomposes the gas velocity into two components 
one of which denotes the Keplerian rotation and the other of which does the 
deviation from it. This scheme enables us to solve the centrifugal balance of a 
gas disk against gravity with better accuracy, since the former inertia force cancels 
the gravity. It is applied to circumbinary disk rotating around binary of which 
primary and secondary has mass ratio, 1.4:0.95. The gravity is reduced artificially 
softened only in small circular regions around the primary and secondary. The 
radii are 7 % of the binary separation and much smaller than those in the previous 
grid based simulations. 7 Models are constructed to study dependence on the gas 
temperature and the initial inner radius of the disk. The gas accretion shows both 
fast and slow time variations while the binary is assumed to have a circular orbit. 
The time variation is due to oscillation of spiral arms in the circumbinary disk. 
The masses of primary and secondary disks increase while oscillating appreciably. 
The mass accretion rate tends to be higher for the primary disk although the 
secondary disk has a higher accretion rate in certain periods. The primary disk 
is perturbed intensely by the impact of gas flow so that the outer part is removed. 
The secondary disk is quiet in most of time on the contrary. Both the primary and 
secondary disks have traveling spiral waves which transfer angular momentum 
within them. 

Subject headings: accretion — hydrodynamics — stars: binary — stars: forma- 
tion 
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INTRODUCTION 



Accretion from a circumbinary disk is expected to play an essential role in binary for- 
mation. Current star formation theory assumes that stars have only 0.01 M Q at their birth 
and gain most of the mass at the main sequence stage through gas accretion. The accretion 
is mainly through disks surrounding the stars. When the stars are members of binaries, they 
are expected to accrete gas from the circumbinary disks when they are young. 

Recent infrared imaging observations support the gas accretion from the circumbi- 
nary disk. The disk ar ound CoKu Tauri/4 is of circumbinary rather than the transitional 
( jlreland fc Kraud 120081 ) . The inner hole of the disk is likely to be cleared by the newly 
discovered binary w ith the separation of ~8 AU. UY Aur has a circumbinary disk having 
a clumpy structure (IHioki et al.l 120071 ). which implies dynamical interaction with the binary 
inside the disk. 

Despite the importance, gas accretion is still mysterious in binar y system. The first nu- 
merica l simulation of the gas accretion onto binary was performed by lArtymowicz fc Lubow 
( 119961 ) . They demo nstrated mass flow throu gh gaps in circumbinary disks by SPH simula- 
tions. Shortly later iBate fc Bonnnelll (119971 ) demonstrated that the secondary has a much 
higher accretion rate than the primary when accreting gas has a relatively large angular 
momentum. Their result seems to be plausible since the secondary has a larger orbit and 
hence is closer to the inner rim of the circumbinary disk. If this is the case, however, the 
mass ratio should increase and approach to unity. Accordingly equal mass binaries should be 
dominant in contradiction to observed sta tistics. Unequal mass binaries are rather common. 
According to iDuauennov fe Mayor! (Il99ll ). the mean mass ratio is ~ 0.4 for the solar type 



star stars. iReid fc Gizij (119971 ) finds that the mass ratio distribution is approximately flat. 
Although the statistic may suffer from observational selection bias, unequal mass binaries 
are surely not rare. Present SPH simulations of star formation predi ct dominance of nearly 
equal mass binaries and do not reproduce the statistics (see, e.g., iDelgad-Doanatei 12004 ; 
Clarke! boosh . 



Ochi et al.l (120051 ). which is referred as Paper I in the following, computed gas accretion 



onto binary using a grid based code and obtained qualitatively different results. The primary 
has a larger accretion rate in their simulations. Their simulations employ different methods 
and assumptions from those of the previous SPH simulations. First the SPH simulations are 
three dimensional, while Paper I is two dimensional. Second the SPH simulations assume 
a cold isothermal gas while Paper I assumes a warm isothermal gas. Third the method 
of gravity softening is different. The SPH particles are absorbed when they enter into 
the regions close to either of the primary and secondary. On the other hand, the gas is 
accumulated in the regions around the primary and secondary in Paper I. The regions of 



- 3- 



artificial gravity softening are much wider in Paper I. The radii are 5 % of the binary 
separation in the SPH simulations while they are 20 % of the binary separation in Paper I. 

The differen ces mentioned above may not account for the difference in the results (see, 



e.g. . IClarkdl2008l ) . The gas is concentrated on the orbital plane in the SPH simulations. Two 



dimensional SPH simulations assu ming a warm g as does not reproduce the gross features 



of Paper I, as shown in Figure 2 of IClarkd (120081 ) . This discrepancy should be resolved to 



clarify the mystery of unequal mass binaries. 

We have developed a new scheme to solve gas accretion onto binary, since an ordinary 
grid-based scheme does not work in practice. 

A standard scheme demands a high spatial resolution around the primary and secondary 
so that the grid spacing is much smaller than the scaleheight. The scale height around the 
primary is evaluated to be Hi = r\c 2 s / '(GM), where r±, c s , G, and Mi denote the radial dis- 
tance from the primary, the sound speed, the gravitational constant, and the primary mass, 
respectively. Then the required spatial resolution is expressed as Ar <C (c a /v v ) 2 r, where 
v v denotes the Keplerian rotation velocity around the primary, a/ G~M\pr\. The required 
resolution is enormous even when the gas is relatively warm. Paper I studies the warm gas 
of c s = 0.25 in the unit system in which the binary separation and angular velocity are 
unity. When the mass ratio is q = 0.7, the scaleheight is Hi = 0.021 at n = 0.2, i.e., at 
the outer edge of the region of artificial gravitational softening. If the gas is colder or closer 
to the primary, the scaleheight is shorter and higher spatial resolution is required. 

We have overcome the difficulty by introducing the reference velocity. The reference 
velocity is set so that the acceleration cancels the gravity. The introduction of reference 
velocity reduces the source term in the hydrodynamical equations. As a result we have 
succeeded in reducing the softening radius down to 7 % of the mean separation. 

The new scheme improved solutions around the primary and secondary. This improve- 
ment is essential to examine the discrepancy between the SPH simulations and Paper I. The 
discrepancy is prominent within the Roche lobe. The accretion from the circumbinary disk 
is mainly through L2 point, i.e., the side close to the secondary. The gas streams around 
the secondary in Paper I while it accretes onto the secondary in SPH simulations. Although 
shock waves form in both types of simulations, the location is scheme dependent. The 
recipient of gas accretion depends critically on the shock location, since the shock energy 
dissipation settles the gas. The flow is highly supersonic in the circumstellar disks, since 
the Keplerian rotation velocity is v v = 1.71 (rx/0.2)- 1 / 2 for q = 0.7. The Mach number 
should exceed 30, if c s = 0.05 as assumed in the SPH simulations. Highly supersonic flow is 
not easy to handle. A small numerical fluctuation may cause spurious turbulence and shock. 
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Furthermore a stationary shock produces a strong density jump since the density enhances 
inversely proportional to the square of the Mach number. Our new scheme resolved spiral 
shock waves formed in the cicumstellar disks. 

Moreover we employed nested grids to obtain higher spatial resolution around the pri- 
mary and secondary. The nested grids consists of rectangular grids having different resolu- 
tions. The highest resolution of 0.12 % binary separation is achieved in the regions close 
to the primary and secondary. On the other hand, the outer boundary is set at 39.2 times 
distant from the center of gravity. The nested grids eliminate possible wave reflection at the 
outer boundary. 

We also modified initial and outer boundary conditions from Paper I. The circumbinary 
disk is assumed to rotate with the Keplerian velocity and to have uniform surface density 
beyond the initial inner edge at the initial stage. The outer boundary condition is set so that 
the surface density and velocity are fixed to be constant in time. These initial and boundary 
conditions reduce fluctuations due to gas accretion. The initial stage may correspond to a 
stage at which accreting gas has been accumulated in the circumbinary disk. 

This paper is organized as follows. Our model and numerical methods are summarized 
in §2. The numerical results are shown in §3 with emphasis on oscillation of spiral waves 
in the circumbinary disk. The oscillation results in time variation in the accretion rate. 
We discuss implications of our results in §4. The technical details on the reference velocity 
method are described in Appendix. 



2. Model 

We consider a binary system in which circumbinary and circumstellar disks are coplanar 
with the orbit. The masses of the disks are neglected and the binary orbit is assumed to 
be circular for simplicity. Furthermore we integrate the hydrodynamical equation in the 
direction vertical to the orbital plane to reduce the problem in the two dimensional. Then 
the hydrodynamical equations are expressed as 

— + V(Su) = 0, (1) 

y + (ii-v)M + -vn = -v$-f]x M (2) 

in the frame corotating with the binary. The symbols, E and II, denote the surface density 
and vertically integrated pressure, respectively. The symbols, u and Q, denote the velocity 
in the corotating frame and angular velocity of the binary, respectively. The center of gravity 
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is located at the origin, r = 0. The gravitational potential (<&) at the position r is expressed 
as 

t x \r-r k \ 2 l 

where and denote the mass and position of the fc-th component star, respectively. For 
simplicity we assume that the gas is isothermal and 

n = c % (4) 

where c s denotes the sound speed. It is set to be either c s = 0.22 or 0.25 in this paper. 

In the following we use the non-dimensional unit system in which the unit length is the 
binary separation, |r*i — r 2 |, and the unit time is the inverse of the angular velocity, 
The gravitational constant, G, is taken to be unit for simplicity. Then the total mass of 
binary is also unity, M\ + M 2 = 1, in this unit system. The orbital plane is assumed to 
coincide with z = in the Cartesian coordinates. The primary is assumed to located at 
ri=(M 2 ,0), while the secondary at r 2 =(— Mi,0). 

We decompose the gas velocity into two components, given (w) and unknown (v), 

u(r,t) = w(r) + v(r,t). (5) 

The former (w) is called reference velocity in the following. Substituting Equation (jSJ) into 
Equation (j2J), we obtain 

dv 1 

— + (u ■ V) v + - VII = g cS -2Qxv-(v V)w, (6) 

where 

g cS = -V$ - (w ■ V) w - 2fl x w. (7) 

The reference velocity is chosen so that the effective gravity, g cS , is small everywhere in 
the computation domain. Thus it is nearly the Keplerian rotation around the each star 
in the vicinity while it is the Keplerian rotation around the center of the gravity in the 
circumbinary disk. The gravitational potential, is softened only in small areas around the 
primary and secondary, i.e., in the regions of \r — Tj| < The softening radius is taken to 
be Pi = P2 = 0.07 in most models. 

Note that the right hand side of Equation does not contain any derivative of v. Thus 
we integrate Equation (jH]) coupled with Equation ([T]) numerically by an ordinary method. 
Although the residual velocity (v) is small, the nonlinear term is taken into account. The 
introduction of the reference velocity is not approximation but purely mathematical trans- 
formation. Further details on the numerical integration are given in the Appendix. 
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The hydrodynamical equations are integrated numerically on the nested grid consisting 
of ten rectangular grids. The rectangular grids contain 512 2 square cells each. Two of 
them have the grid spacing of Ax = 1.2 x 10~ 3 and cover the vicinities of the primary 
and secondary. Another rectangular grid has the grid spacing of Ax = 0.1536 and covers 
the whole computation region of 78. 6 2 square area. The rest grids have intermediate spatial 
resolutions and cover a part of the computation region. They are arranged so that the spatial 
resolution is higher in the regions closer to the primary and secondary. 

The outer boundary is placed at r = 38.9. The surface density and velocity are replaced 
with the initial values at each time step outside the outer boundary. The outer boundary is 
so far from the binary that it has no appreciable effects on the simulations. 

The initial surface density is taken to be 

= E out + S in E out ~ S in / r^ro\ 

2 2 V h J 

where S out = 1.0, S in = 0.1, and h = 0.144. Thus the initial circumbinary disk has an 
inner edge at r = r . The initial velocity is set so that the residual velocity, v, vanishes. 

3. RESULTS 

This paper shows 6 models having different c s , r , and (3\ (= $2). The model parameters 
are listed in Table [TJ First we introduce model 1 because it has the smallest c s and a modestly 
large Tq. When c s is smaller (i.e., the gas is colder), the disk is geometrically thinner and our 
2D approximation is better. When ro is larger, the cirucmbinary disk has a larger specific 
angular momentum and the gas can accrete onto the binary only after substantial angular 
momentum transfer. 



3.1. Model 1 

The circumbinary disk has the inner edge at r = 2.5 at the initial stage of model 1. 
Figure [1] shows the surface density distribution at t — 15, 30, 90, and 288 in model 1. The 
color denotes the surface density in logarithmic scale while the arrows denote the velocity in 
the frame corotating with binary, u. The primary (right) has a more massive circumstellar 
disk than the secondary (left). 

The circumstellar disks accretes gas mainly through L2 point from the circumbinary 
disk. This is reasonable since the Roche potential is lower at L2 point by A$ = 6.70 x 10 -2 . 
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The potential difference is appreciably larger than the square of the sound speed, c 2 s = 
4.84 x 10 -2 . The circumbinary disk has two-armed spiral shock waves (denoted by yellow in 
Fig.d]). 

The gas velocity is small near L2 point. We surmise that it is regulated to be the sound 
speed, since the Roche potential has a saddle point at L2. Then L2 point should be similar 
to the throat of Laval nozzle. As well known, the gas velocity coincides with the sound speed 
at the throat. Similarly the gas velocity should be regulated at the sound speed at L2 point. 

Figure [2] is a close up view of model 1 at t — 90.0. The left panel is focused on the 
regions of high surface density while the right is on the regions of low surface density. Spiral 
shock waves are excited also in the circumstellar disks. The outer part of the primary disk is 
disturbed appreciably while the secondary is not. The disturbance is likely due to the impact 
of gas flow from L2. The gas stream is nearly normal to to the primary disk while it is nearly 
tangential to the secondary disk. See Figure [3] for the radial velocity toward each star. The 
primary disk receives a larger impact, which produces 'a hot spot' on the side closer to LI 
point. The radial gas inflow is shown as a function of the azimuthal angle centered on the 
primary in Figure 0} The ordinate denotes the mass flux, i.e., the product of the surface 
density and the radial velocity toward the primary at t — 90.0 and |t* — r*i| = 0.23. The 
azimuthal angle is (j) — and ir/2 in the directions of increasing x and y, respectively. 

Figure shows the region of large V • u at t — 90.0 by color. The colored is the shock 
front, since the gas is compressed strongly there. It is clearly shown that the primary disk 
is associated with the strong shock while the secondary disk is not. Furthermore the gas 
stream changes its density and velocity semi-periodically on a short timescale. The impact 
on the primary disk varies accordingly. 

Note that Figure [2] is d ifferent from both the cold and warm SPH models of Delgado- 



Donate shown in Figure 2 of lClarkd (120081 ). The secondary is surrounded by a shock wave in 
his cold SPH model while no shock wave is seen in his warm SPH model. The shock location 
is closely related to the mass accretion rate of each component star. It is crucial whether 
the gas streaming from LI point causes a shock around the secondary. 

The primary disk mass increases while oscillating appreciably. Figure [6] shows the 
primary disk mass as a function of time. Each curve denotes the total mass of gas enclosed in 
a circle around the primary. It also indicates that the surface density falls off sharply around 
\r — Ti| = 0.2. Thus we regard the total mass enclosed in the circle of \r — rx\ = 0.2 as 
the disk mass in the following. 

Figure [7] is the same as Figure [6] but for the secondary disk. The secondary disk mass 
is smaller than the primary one and increases steadily. The gas flow around the secondary 
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is smooth and the secondary disk accretes gas without significant disturbance. 

Semi-periodic change is also seen in the circumbinary disk. Figure [H] shows the time 
variation in the surface density measured at four points on the axis of y = 0. The period of 
surface density variation is about P ~ 4.4 when measured at a given point in the corotation 
frame. Variation is also seen in the radial velocity. The variation has the temporal and 
azimuthal dependence of cos [2(0 — ut)] in the corotation frame, where the angular frequency 
is approximately uo ~ 0.7f2 in the period 20 < t < 80. Around t ~ 100, the surface density 
variation is not sinusoidal; two waves having similar but different angular frequencies seem 
to be excited. The amplitude of the wave varies on a longer timescale as shown in Figure [HJ 
It is large in the period t > 130. This large variation is due to one armed low surface density 
region shown in the lower right panel of Figure [TJ 

This quasi-periodic change has not been reported in the literature. It may be missed or 
smeared out in the previous simulations by some reasons. Paper I and some other simulations 
assumed continuous dynamical accretion ont o the circumbina r y disk , whic h grows the disk 



and m ay hinder our perception. Although iGiinther fc Kleyl (120021 ) and iGiinther fc Kiev 



(120041 ) assumed quasi-static initial and boundary conditions, they employed a large empirical 
viscosity and spatial resolution much lower than ours. The numerical viscosity and relatively 
low spatial resolution may damp down the quasi-periodi c change. We have learned that 



recent SPH simulations have shown similar variation from I Clarke! (120091 ) 



Figure [2] shows a thin stripe of high surface density connecting primary and secondary 
disks as in Paper I. The stripe is formed by the collision of the gas flows rotating around 
the primary and secondary. In practice, it is the boundary between the domains in one of 
which the gas flows from secondary to primary and in the other of which it flows vice versa. 
It is located around LI point but its location varies with time. When the inflow from L2 
has a higher surface density, the stripe is closer to the primary. The location of the stripe is 
critical to which component of the binary accretes more. 

The sum of the disk masses amounts to Mid + A^2d = 1-78 at t = 144. This is equal to 
the mass of the circumbinary ring of Ar = 0.11 at r = r . The total disk mass amounts to 
Mi d + M 2d = 4.75 at t = 288. 

The time derivative of a disk mass gives us the accretion rate in principle. It is, however, 
highly variable and hard to read gradual change in the accretion rate. To erase out the short 
time variability, we obtain the average accretion rate for time intervals of At = 2.88. Figure^ 
shows the average accretion rates for the primary (solid), dM ld /dt and secondary (dashed), 
dM2d/dt. The accretion rate is very high in the first a few rotation period (t < 25). The 
primary has a higher accretion rate than the secondary except around t ~ 230 and 270. 



- 9- 



It should be noted that there is correlation between the accretion rates and oscillation in 
the circumbinary disk. The accretion rates are low when the circumbinary disk is relatively 
quiet, i.e., around t ~ 100. This indicates that the oscillation enhances the accretion from 
the circumbinary disk through angular momentum exchange. 

Figure [10] shows the azimuthally averaged surface density in the circumbinary disk as 
a function of r and t. It shows repeatedly excited waves propagating outwards of which 
phase velocity is dr/dt ~ 0.25 and 0.24 at r = 3.2 and 4.2, respectively. The inner edge of 
the circumbinary disk retreats very slowly. Figure [TD] shows that the amplitude of the wave 
varies on the timescale of several tens rotation period. 

At some epochs, the secondary has a higher accretion rate than the primary. Figure [TT1 
shows the flow at t — 230.4 as an example of such a stage. As shown in the figure, the outer 
part of the circuprimary disk is shed to flow into the secondary lobe. The circumsecondary 
disk receives the flow and has a hot spot on the LI point side. Similar mass transfer is also 
seen around t ~ 150 and 190. The inflow is still mainly through L2 point and the inflow 
through L3 point is minor. Note that the gas flows from primary to secondary around LI 
point in this period. Shortly after this stage, the circumsecondary disk is disturbed very 
much and captures a gas flow from L2 point. The gas flow inside the Roche lobe is crucial 
for the branching of the inflow from the circumbinary disk. 

3.2. Model 2 

We compare model 2 with model 1 to examine the dependence on r , the initial inner 
radius of the cicumbinary disk. The model parameters of models 2 are the same as those of 
model 1 except ro- 

Figure [12] shows the surface density distributions at t — 162. The circumstellar disk 
is more massive in model 2 since the circumbinary disk has a smaller inner radius. The 
primary disk is more massive than the secondary one also in model 2. The primary disk 
mass increases while oscillating appreciably as shown in Figure [TBI The secondary disk mass 
increases slowly and smoothly as shown in Figure HH The averaged accretion rate is shown 
in Figure [T5l 

3.3. Models 3, 4, and 5 

We compare model 4 with model 1 to examine the dependence on the gas temperature. 
The gas temperature is 29.1% higher in model 3 than in model 1 since the isothermal sound 
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speed is set to be c s = 0.25. It is still lower than the the potential gap between L2 and L3 
points. Thus the flow is relatively cold although it is quite higher than that assumed in SPH 
simulations. 

Figure [TBI shows the distributions of £ and u at t = 108 and 216 in model 3. The spiral 
arms have a little wider opening angle in the circumbinary disk of model 3 than in that of 
model 1. The gas inflows mainly from L2 point and the surface density is low around L3 
point. Remember that the potential difference is A$ = 6.70 x 10 -2 between L2 and L3. It 
is larger than c 2 = 6.25 x 10 -2 also in models 3, 4, and 5. 

Figure [T71 compares the disk masses of model 3 with those of model 1. Both the primary 
and secondary disk masses are ~50 % larger in model 4. This is an envidence that the 
oscillating spiral arms are responsible for the accretion from the circumbinary disk. When 
c s is higher, the angular momentum exchange should be more effective. 

Figure [TBI is the same as Figure [TTl but for comparison between models 2 and 5 in which 
only the sound speed, c s , is different. The disk masses are ~30 % larger in model 5. When 
the gas temperature is higher, the accretion rate is higher. The gas temperature affects a 
little the ratio between the primary and secondary accretion rates. 

The vibration of the spiral shock wave is also seen in model 3. Figure [19] is the same as 
Figure [8] but for model 3. The oscillation frequency is almost the same as that in model 1. 

3.4. Model 6 

We have constructed model 6 to examine the effect of the softening radius. The initial 
state is the same as that of model 3 except for f3\ and (3%. The gravity is artificially reduced 
in the regions of \r — ri\ < 0.14 and \r — r%\ < 0.14 in model 6. 

Figure I2"01 shows the surface density distribution at t = 90.0 in model 6. Also in model 6, 
the the primary disk is more massive than the secondary one. Since the gravity is artificially 
reduced in a larger area, the primary and secondary disks are rings with larger holes in the 
early stages. It takes a longer time for the gas rings to spread into the center. 

Figures [UJ and [221 are the same as Figures [6] and [7] but for the primary and secondary 
disks in model 6, respectively. The mass of the secondary disk is smaller than in model 3. 
Note, however, that the softening radius is still as small as /?2 = 0.14. The secondary disk 
extends outside the softening radius and can capture gas if it is accreted. The primary disk 
sheds its outer part at some epochs (e.g. t ~ 150 and 190) and the secondary receives it. 
This feature is also seen in model 6. 
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The softening radius affects the inner structure of the circumstellar disks. However, it 
does not affect the gas stream from L 2 point directly; the gas flows far from the region of 
artificial reduction in the gravity. We conclude that the softening radius does not give a 
serious effect on the result that the primary accretes more than the secondary. 

3.5. Model 7 

We have constructed model 7 to examine the effect of initial surface density distribution. 
Model 7 is the same model 2 except for h which denotes the surface density gradient at the 
initial disk edge. The surface density increases gradually in model 7 since h = 0.6. Hence 
the initial pressure force is much weaker than in model 2, since h = 0.144 in model 2. 

Figures [23] and [24] show the disk mass as a function of time in model 7 for primary and 
secondy, respectively. The initial rise in the primary disk mass is steeper in model 7. This 
is because the initial disk is extended more in model 7. Inner part of the circumbirnary disk 
is truncated even when the pressure gradient is weak. The primary disk mass of model 7 is 
almost the same as that of model 2. The disk mass does not depend seriously on h. 

4. DISCUSSIONS 

In this section, we discuss the softening radius, the hot spot, the accretion impact on 
the circumstellar disk, and the oscillation of spiral waves excited in the circumbinary disk. 

In this work the softening radii, f3\ and {3%, are three times smaller than those in our 
previous work (Paper I). They are also much smaller than the circumstellar disk radii. Our 
numerical simulations clearly resolved internal structures of the circumstellar disks such as 
spiral waves within them. The simulations still depend a little on the softening radius. The 
gas accretion is slow inside the softening radii and the circumstellar disks have inner holes in 
the early period. The circumstellar disks affect accretion rates through dynamical interaction 
with the gas stream from L2 point. 

Further reduction in the softening radius requires extremely high spatial resolution and 
accordingly unfeasible computation cost. See Appendix on the spatial resolution required. 

As shown in the previous section, the primary disk has a hot spot during the primary 
accretion rate dominates over the secondary one. The hot spot is due to the impact of the gas 
stream and the gas looses its kinetic energy through shock dissipation. The kinetic energy 
dissipation plays an important role for a gas element to be accreted onto a circumstellar 
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disk. 

It is well known that the Jacobi integral, 



\u\ 2 

J=^ + $, (9) 

is a constant of motion for a particle moving around a circular binary. After some algebraic 
manipulation, we obtain 

DJ' 9 «91n£ , . 

-m = (10) 

for an isothermal gas flow without shock, where 

J 1 = J + c \ InS. (11) 

The difference between J and J' is small when the gas temperature is low. The right hand 
side of Equation fflOl) is also small. Thus the Jacobi integral can be regarded as a constant 
for a Lagrangian gas element if the gas is cold. When a particle orbits around the primary 
at the Keplerian velocity, the Jacobi integral is evaluated to be 




Jl = _ Ml __ 5i + 5^ F= ^ r |r-r 1 |j , (12) 

where the tidal force due to the secondary is neglected for simplicity. Similarly, it is evaluated 
to be 

+K\&'H 2 ' (i3) 

for a particle orbiting around the secondary. The values of J\ and J 2 are shown as a function 
of the orbital radius in Figure The dashed lines denote the potential energy at the 
Lagrangian points L 1; L 2 , and L 3 . 

The inflowing gas has a small velocity in the corotation frame when it passes L 2 point. 
Thus the Jacobi integral is nearly equal to the potential energy at L 2 point. Provided that the 
gas accretes onto either the primary or the secondary without changing the Jacobi integral. 
Then the orbital radius should be as much as \r — r\\ = 0.34 or \r — r 2 | = 0.30. The orbital 
radius can be smaller only when the Jacobi integral is reduced substantially through strong 
shock dissipation. 

When the orbital radius is \r — t*i| = 0.26, the Jacobi integral, J\, is equal to that at rest 
at Li point. This implies that a gas element is marginally bound to the primary when the 
orbital radius is 0.26. The gas element may flows into the secondary lobe if perturbed a little. 
Our simulations indeed show such mass transfer from the primary disk to the secondary one 
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around t ~ 190 in model 1. The outer part of primary disk is shed into the secondary. 
The radius of the marginally bound orbit is \r — r 2 \ = 0.21 for the secondary disk. The 
secondary disk is confined mostly inside the marginal orbit in our simulations. 

The above discussion gives a guideline for choosing the softening radius. A gas element 
orbiting at the softening radius should have an appreciably smaller Jacobi integral than 
resting at LI point. If the gas element has an eccentric orbit around a star, the pericenter 
should be located outside the softening radius. Since the softening radius was 20 % of the 
binary separation in Paper I, it is close to the radius of the marginally bound orbit for the 
secondary and should be set smaller. 

The kinetic energy dissipation may brighten the hot spot at some wavelengths. Although 
the hot spot is assumed to be isothermal in our simulations for simplicity, it should be hotter 
than the rest part of the circumstellar disk. It s hould be ana logous to the hot spot in the 



cataclysmic variables and X-ray binaries (see, e.g.. lMarshll200ll ). In these systems, a compact 
star (either of white dwarf, neutron star and black hole) accretes gas from its companion 
through Li point. The gas inflow forms a hot spot on the circumstellar disk around the 
compact star. It is often identified by occultation. Similarly, the hot spot may be observed 
in binaries of young stellar objects. It may be identified by the Doppler effect if it is luminous 
at some emission lines. 

It should be noted that the impact depends on the circumprimary disk. When the 
circumprimary disk has a higher surface density and pressure, the inflowing gas stream is 
relatively weaker and hence produces the impact point at an outer radius. The evolution of 
circumstellar disks is controlled by angular momentum transfer and hence will be affected 
by assumed viscosity and gravity softening. 

Next we discuss the oscillation of spiral waves excited in the circumbinary disk. The 
oscillation is closely related to gas accretion from the circumbinary disk as shown in the 
previous section. This is also confirmed by Equation (fTO]) . If the density distribution is 
stationary in the corotation frame, the modified Jacobi constant, J', is constant except 
when it gets across a shock front. The shock dissipates the kinetic energy and decreases 
the Jacobi constant. If the Jacobi constant decreases, the orbital radius increases. This is 
because the Jacobi constant is evaluated to be 

J =-Tr~^' (14) 

when the gas element rotates around the binary with the orbital radius, r, in the circumbinary 
disk. Equation (TH|) is derived under the assumption that the orbit is almost circular and 
the velocity is Keplerian. The Jacobi constant is larger for a larger r when r > 1. In other 
words, the gas element gains angular momentum from the spiral wave corotating with the 
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binary, since the angular velocity of the binary is higher than that of a gas element in the 
circumbinary disk. This means that a stationary shock wave drives the circumbinary disk 
outward. The accretion is due to vibration of the shock wave. Some gas elements move 
inward by loosing the angular momentum and Jacobi constants through exchange with some 
other gas elements, when the spiral shock wave is not stationary. 

The oscillation has a typical period of ~ 4.5 when measured in the corotation frame. 
Accordingly the wave has an apparent phase velocity of ~ 0.7 in the corotation frame since 
the azimuthal wavenumber is m = 2. Then the phase velocity is ~ 0.3 in the rest frame 
and coincides with the angular velocity of gas orbiting around the binary at r ~ 2.2. This 
implies that the wave is excited by resonance at the inner edge of the circumbinary disk. 

The amplitude of the oscillation varies on the timescale of several tens rotation period. 
The variation is coherent in the radial direction. When the amplitude is small, the whole 
cicumbinary disk is nearly stationary. When the amplitude is large, both the radial and 
azimuthal components of velocity oscillate at the same frequency. The variation may induce 
variability in the accretion onto the binary. 

We thank Tomoaki Matsumoto for his useful comments and advices. This study is 
financially supported in part by the Grant-in-Aid for Scientific Research on Priority Area 
(19015003) of The Ministry of Education, Culture, Sports, Science, and Technology (MEXT). 



We describe the reference velocity methods for integrating hydrodynamical equations. 

First we divide the computation domain into four areas: one around the primary, another 
around the secondary, another surrounding the binary, and the rest. The reference velocity 
is defined differently in each area. 

The reference velocity is defined to be the Keplerian rotation around the center of mass, 
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REFERENCE VELOCITY METHOD 



w = oo(\r\) e z x r 



(Al) 




1 + r" 3 / 2 , 



(A2) 



outside the circle of \r\ > 1. 



The reference velocity is taken to be circular rotation around the primary, 




(A3) 
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in the vicinity of \r — r\\ < 71. The angular velocity is set to be 



1/2 



u(\r - t*i|) = < 



n + 



Mi 



7 3 / \r — Ti| 

4 _ 4 

1/2 
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r — t*i| 

' 5 3 |r — 7*1 1 
2 ~ 2 



Mi 



(k - ri| < ft) 

(ft < |r - rij < «i) 

(ai < |r - ri| < 71) 



(A4) 

The angular velocity coincides with the Keplerian velocity around the primary in the region 
of /3i < ri < ai. It is slower than the Keplerian velocity in the region of a\ < T\ < 7i- 
The radius, 71, is set to be 



Ti 



2ai 



(A5) 



so that the angular velocity vanishes at \r — r±\ = 71. 



The gravity of the primary is softened so that it balances with the centrifugal force due 
to the reference velocity. Thus the gravity due to the primary is reduced to 



g (\r - n|) = --^-3 (r - r 



7 |r — 7*1 1 
4 



ft 2 



(A6) 



Similarly the reference velocity is taken to be circular rotation around the secondary in 
the vicinity of the secondary. It is expressed as 



w = u(\r - r 2 \)e z x (r - r 2 ) 



(A7) 
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4 ~ 4 
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where 
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a: 
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3 |r - r 2 | 
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(|r - r 2 | < ft) 

(ft < |r - r 2 \ < a 2 ) 
(a 2 < \r — r 2 \ < 72) 



2a 2 



5 /Mo 



2 V«l 







(A? 



(A9) 



The values of a iy ft, and 7$ are tabulated in Table [2j Note 71 + 7 2 < 1 and accordingly 
that the region of |r — 7*1 1 < 71 does not overlap with that of \r — r 2 \ < 7 2 . 
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The reference velocity vanishes in the rest area, i.e., inside the circumbinary disk but 
outside the primary and secondary disks. Remember that the reference velocity is continuous 
even on the boundary between the different areas. Thus the derivative of the reference 
velocity is regular and continuous except discontinuity on the boundary between the areas. 

The discontinuity of the reference velocity derivative has no serious effect on the evalu- 
ation of the source term. The reference velocity is analogous to the vector potential in the 
electromagnetism. Thus its derivative is analogous to the magnetic field. The magnetic field 
can be discontinuous so the derivative of the the reference velocity can be. Remember that 
the reference velocity is divergence free, V • w = 0. This is analogous to the Coulomb gauge 
in the electromagnetism. 

The numerical flux is evaluated at the cell surface as in the ordinary finite volume 
method. The reference velocity is evaluated analytically at the cell surface while the other 
variables are evaluated from the cell center values by using interpolation with the minmod 
limiter. Then the reference velocity is numerical viscosity free and the numerical solution 
has a larger Reynolds number. 

Introduction of the reference velocity improves evaluation of the centrifugal force greatly. 
The centrifugal force, v • Vti, contains a spatial derivative of the velocity. The spatial 
derivative is evaluated from the spatial difference in an ordinary finite difference method 
and truncation error is inevitable in its numerical evaluation. The truncation error should 
be small so that minor forces such as pressure force are appreciated properly. Note th at the 



nume rical derivative is only the first order accurate near shock fronts in a TVD scheme (IToro 



20091 . see, e.g.,). When and only when the spatial resolution, Ax, is much smaller than the 



the scaleheight, 

A* « ^ (A10) 

the truncation error is ensured to be much smaller than the pressure force. This condition 
is equivalent to 

Ax < ( — J r , (All) 



vk 



where vk denotes the Keplerian velocity. 



The reference velocity relaxes the above mentioned requirement on the truncation error. 
Alternatively, the resolution should be so small that the difference in reference velocity is 
smaller than the sound speed between any adjacent cells. Otherwise a gas element would be 
accelerated (or decelerated) too much by advection from a cell to another. This requirement 
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is less tight, since it can be expressed as 

Ax « r . (A12) 

Equation flAllj) is much more tight than this, when the Mach number (vk/c s ) is large. The 
Mach number is as high as 

/ S'f / | |\ -1/2 

V K 100/ M\ \ / r — ri' x 



Ifl = 13.3 ^ LiL (A13) 

c s V0.596 y V 0.07 J V0.22/ v ' 

for the Keplerian rotation around the primary. Thus the reference velocity. Further reduction 
in the softening radius requires extremely high spatial resolution even when the reference 
velocity is introduced. 

Our numerical scheme is sim ilar to FARGO (Fast Advection in Rotating Gaseous Ob- 



jects) developed by iMassetl (120001 ) for solving hydrodynamics a gaseous disk rotating around 
a star on fixed polar grids. FARGO decomposes the rotation velocity into two components, 
the azimuthally averaged rotation and deviation from it. The advection due to the former 
is treated separately in FARGO. The separation relaxes the ordinary Courant condition and 
speeds up the computation. At the same time, numerical diffusivity is much smaller in 
FARGO thanks to the two step advection. The centrifugal force is dominated by that due 
to the averaged rotation, J 2 /r 3 , where J denotes the specific angular momentum. Thus the 
truncation error is reduced also in FARGO. 

The reference velocity is, however, different from FARGO. Application of FARGO is 
limited to the case that the grid is set almost parallel to the gas flow. The reference velocity 
method can be applied to any numerical grids. FARGO can take a long time step while the 
reference velocity method can not. 



Our numerical scheme is also similar to that of iLe Vequd (119981 ). He proposed a scheme 
to take account of gradient of a physical variable within a numerical cell for elliminating 
source terms apparently. In case of hydrodynamics, the gravity can be taken into account 
if the pressure has different values at the cell surfaces opposing each other. Our reference 
velocity has different values at the cell surfaces to incorporate the gravity. His scheme is 
effective for a flow in a quasi-static balance while ours is effective for a flow supported mainly 
by the centrifugal force. 
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Fig. 1. — Each panel shows the surface density (color) and velocity distributions at t — 15, 
30, 120, and 240 in model 1. The dashed curves are the contours of the Roche potential. 
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Fig. 2. — The surface density (color) and velocity (arrows) distributions inside the Roche 
lobe at t — 90.0 in model 1. The dashed curves are the contours of the Roche potential. The 
left and right panels denote the same stage but in the different color scales. 
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Fig. 3. — The color denotes the radial velocity toward each star at t — 90.0 in model 1. 
The left panel denotes the velocity toward the primary, while the right does that toward the 
secondary. 
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Fig. 4. — The solid curve denotes the radial mass flux, £i> r i, onto the primary disk as 
a function of the azimuthal angle, <f>i = tan _1 [t//(:r — x±)], where v r \ denotes the radial 
velocity with respect to the primary. The dashed curve is the same as the solid one but for 
the secondary. 
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Fig. 5. — The region of strong shock compression is marked by color on the surface density 
map at t = 90 in model 1. The colored region is identified by the amplitude of V ■ u. It is 
considerably broadened artificially, since the the shock wave is very sharply resolved and is 
hard to be recognized in the raw data. 
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Fig. 6. — Each curve denotes the total mass of gas inside a circle around the primary as a 
function of time. The numbers denote the radii of the circles. 
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Fig. 8. — The time variation of the surface density on the axis of y = in model 1. 
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Fig. 9. — The accretion rates of the primary and secondary derived from the best fit poly- 
nomials to their disk masses. 
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Fig. 10. — The evolution of the azimuthally averaged surface density of the circumbinary 
disk in model 1. 
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Fig. 11. — The same as the left panel of Fig. [2] but for the stage of t = 230.4. 
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Fig. 12. — The surface density distribution around the binary at t — 162.0 in model 2. 



31 




0.04 
0.03 

£ 0.02 

-o 0.01 

0.00 
-0.01 



Primary 




Jl 



Secondary 




ji 



50 100 150 200 250 
t 



Fig. 15. — The same ad Fig. [9] but model 2. 
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Fig. 16. — The surface density (color) and velocity (arrows) distributions at t = 108 and 216 
in model 3. 
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Fig. 17. — The comparison of disk masses between model 1 (thin) and 3 (thick). The solid 
curves denote the primary disk mass while the dashed do the secondary one. 
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Fig. 19- 



The same as Fig. [8] but for model 4. 
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Fig. 20. — The surface density distribution at t — 90.0 in model 6. The notation is the same 
as that of Fig. [TJ 
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Fig. 25. — The Jacobi integral is shown for circular orbits around the primary (solid) and 
secondary (dashed). The abscissa is the orbital radius. The value of Jacobi integral is also 
shwon for particles at rest on the Lagrangian points, Li, L 2 , and L 3 for reference. 
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No. 


C s 




H 


Pi (= &) 


1 


0.22 


2.5 


0.144 


0.07 


2 


0.22 


2.1 


0.144 


0.07 


3 


0.25 


2.5 


0.144 


0.07 


4 


0.25 


2.3 


0.144 


0.07 


5 


0.25 


2.1 


0.144 


0.07 


6 


0.22 


2.5 


0.144 


0.14 


7 


0.22 


2.1 


0.600 


0.07 



Table 1: Summary of models. 



i 


1 


2 


Mi 


0.5957 


0.4043 


OLi 


0.4000 


0.2700 


li 


0.5793 


0.41028 



Table 2: Model parameters for the reference velocity. 



